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Abstract 

The Landau-Lifshitz Navier-Stokes (LLNS) equations incorporate thermal fluctua- 
tions into macroscopic hydrodynamics by using stochastic fluxes. This paper examines 
explicit Eulerian discretizations of the full LLNS equations. Several CFD approaches 
are considered (including MacCormack's two-step Lax-Wendroff scheme and the Piece- 
wise Parabolic Method) and are found to give good results (about 10% error) for the 
variances of momentum and energy fluctuations. However, neither of these schemes 
accurately reproduces the density fluctuations. We introduce a conservative centered 
scheme with a third-order Runge-Kutta temporal integrator that does accurately pro- 
duce density fluctuations. A variety of numerical tests, including the random walk 
of a standing shock wave, are considered and results from the stochastic LLNS PDE 
solver are compared with theory, when available, and with molecular simulations using 
a Direct Simulation Monte Carlo (DSMC) algorithm. 
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1 Introduction 



Thermal fluctuations have long been a central topic of statistical mechanics, dating back 
to the light scattering predictions of Rayleigh (i.e., why the sky is blue) and the theory of 
Brownian motion of Einstein and Smoluchowski [I]. More recently, the study of fluctuations 
is an important topic in fluid mechanics due to the current interest in nanoscale flows, with 
applications ranging from micro-engineering (2j [3l H] to molecular biology (5j El Ej. 

Microscopic fluctuations constantly drive a fluid from its mean state, making it pos- 
sible to probe the transport properties by fluctuation-dissipation. This is the basis for light 
scattering in physical experiments and Green-Kubo analysis in molecular simulations. Fluc- 
tuations are dynamically important for fluids undergoing phase transitions, nucleation, hy- 
drodynamic instabilities, combustive ignition, etc., since the nonlinearities can exponentially 
amplify the effect of the fluctuations. 

In molecular biology, the importance of fluctuations can be appreciated by noting 
that a typical molecular motor protein consumes ATP at a power of roughly 10~ 16 watts 
while operating in a background of 10~ 8 watts of thermal noise power, which is likened to 
be "as difficult as walking in a hurricane is for us" [6]. While the randomizing property 
of fluctuations would seem to be unfavorable for the self-organization of living organisms, 
Nature has found a way to exploit these fluctuations at the molecular level. The second 
law of thermodynamics does not allow motor proteins to extract work from equilibrium 
fluctuations, yet the thermal noise actually assists the directed motion of the protein by 
providing the mechanism for overcoming potential barriers. 

Following Nature's example, there is interest in the fabrication of nano-scale devices 
powered by [8] or constructed using [9] so-called "Brownian motors." Another application 
is in micro-total-analytical systems (//IAS) or "lab-on-a-chip" systems that promise single- 
molecule detection and analysis pjj] . Specifically, the Brownian ratchet mechanism has been 
demonstrated to be useful for biomolecular separation [Til EG] and simple mechanisms for 
creating heat engines driven by non-equilibrium fluctuations have been proposed [131 EH]- 
Finally, exothermic reactions, such as in combustion and explosive detonation, can depend 
strongly on the nature of thermal fluctuations [151 EE] • 

To incorporate thermal fluctuations into macroscopic hydrodynamics, Landau and 
Lifshitz introduced an extended form of the Navier-Stokes equations by adding stochastic 
flux terms [17]. The Landau-Lifshitz Navier-Stokes (LLNS) equations may be written as 

U t + V-F = V-D + V-S (1) 
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where v is the fluid velocity, P is the pressure, T is the temperature, and r = r^(Vv + 
Vv T — |IV ■ v) is the stress tensor. Here rj and k are coefficients of viscosity and thermal 
conductivity, respectively, where we have assumed the bulk viscosity is zero. 

The mass flux is microscopically exact but the other two flux components are not; 
for example, at molecular scales heat may spontaneously flow from cold to hot, in violation 
of the macroscopic Fourier law. To account for such spontaneous fluctuations, the LLNS 
equations include a stochastic flux 
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where the stochastic stress tensor S and heat flux Q have zero mean and covariances given 
by 

(Sy(r, 0) = 2^T {5? k 6« + 5«5? k - §tfg*g) S(r - v')S(t - t% (6) 



and 



(Qi(r,t)Q i (r',0> = 2k B KT 2 5«5(r - r')5(t - t') 



where ks is Boltzmann's constant. The LLNS equations have been derived by a variety of 
approaches (see [T71 [181 EH 120]) and have even been extended to relativistic hydrodynam- 
ics [21]. While they were originally developed for equilibrium fluctuations (see Appendix A), 
specifically the Rayleigh and Brillouin spectral lines in light scattering, the validity of the 
LLNS equations for non-equilibrium systems has been derived [22] and verified in molecular 
simulations I231ET 
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In this paper we investigate a variety of numerical schemes for solving the LLNS equa- 
tions. For simplicity, we restrict our attention to one- dimensional systems, so (p]) simplifies 
to 
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with a being the surface area of the system in the yz-plane. 

Furthermore, we take the fluid to be a dilute gas with equation of state P = pRT 
and energy density E = c v pT + ^pu 2 . The transport coefficients are only functions of 
temperature; for example, for a hard sphere gas rj = t]q\/T and k = Kq\/T, where r] and k 
are constants. The numerical schemes developed in this paper may readily be formulated for 
other fluids. Our choice is motivated by a desire to compare with molecular simulations (see 
Appendix B) of a monatomic, hard sphere gas (for which R = ks/m and c v = -^j- where m 
is the mass of a particle and the ratio of specific heats is 7 = |). 

Several numerical approaches for the Landau-Lifshitz Navier-Stokes (LLNS) equa- 
tions, and related stochastic hydrodynamic equations, have been proposed. The most suc- 
cessful is a stochastic lattice-Boltzmann model developed by Ladd for simulating solid-fluid 
suspensions [25J. This approach for modeling the Brownian motion of particles was adopted 
by Sharma and Patankar [26] using a finite difference scheme that incorporates a stochastic 
momentum flux into the incompressible Navier-Stokes equations. By including the stochastic 
stress tensor of the LLNS equations into the lubrication equations Moseler and Landman [27] 
obtain good agreement with their molecular dynamics simulation in modeling the breakup 
of nanojets. An alternative mesoscopic approach to computational fluid dynamics, based 
on a stochastic description defined by a discrete master equation, is proposed by Breuer 
and Petruccione [28l [29]. They show that the structure of the resulting system recovers the 
fluctuations of LLNS. 
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Serrano and Espanol [30] describe a finite volume Lagrangian discretization of the 
continuum equations of hydrodynamics using Voronoi tessellation. Casting their model into 
the GENERIC structure [31] allows for the introduction of thermal fluctuations yielding a 
consistent discrete model for Lagrangian fluctuating hydrodynamics. Fabritiis et al. [32J S3] 
derive a similar mesoscopic, Voronoi-based algorithm using the dissipative particle dynamics 
(DPD) method. The dissipative particles follow the dynamics of extended objects subject 
to hydrodynamic forces, with stresses and heat fluxes given by the LLNS equations. 

In earlier work Garcia, et al. [31] developed a simple finite difference scheme for the 
linearized LLNS equations. Though successful, that scheme was custom-designed to solve a 
specific problem; it cannot be extended readily, since it relies on special assumptions of zero 
net flow and constant heat flux and would be unstable in the more general case. Related 
finite difference schemes have been demonstrated for the diffusion equation [35], the "train" 
model [35], and the stochastic Burgers' equation [37], specifically in the context of Adaptive 
Mesh and Algorithm Refinement hybrids that couple particle and continuum algorithms. 

In the next section we develop three stochastic PDE schemes based on standard CFD 
schemes for compressible flow. The schemes are tested in a variety of scenarios in sections 
[3] and HI measuring spatial and time correlations at equilibrium and away from equilibrium. 
Results are compared to theoretically derived values, and also to results from DSMC particle 
simulations (see Appendix B). We also examine the influence of fluctuations on shock drift, 
comparing results from the LLNS solver with DSMC simulations. The concluding section 
summarizes the results and discusses future work, with an emphasis on the issues related to 
using the resulting methodology as the foundation for a hybrid algorithm. 

2 Numerical Methods 

The goal here is to develop an Eulerian discretization of the full LLNS equations, representing 
an extension of the approach discussed in [37] to compressible flow. We restrict consideration 
here to finite- volume schemes in which all of the variables are collocated, so that the resulting 
method can form the basis of a hybrid method in which a particle description (DSMC) 
is coupled to the LLNS discretization. Within this class of discretizations, our aim is to 
recover the correct fluctuating statistics. In this section we develop two methods based on 
CFD schemes that are commonly used for the Navier-Stokes equations. We then introduce 
a specialized centered scheme designed to capture fluctuation intensities. 
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2.1 MacCormack Scheme 

Based on the success of the simple second-order scheme in [34] . we first consider MacCor- 
mack's variant of two-step Lax-Wendroff for solving fluctuating LLNsJll The MacCormack 
method is applied in the following way: 

= U? - £ (B? - F?-i) + ^ (D? +1/2 - D?_ 1/2 ) 

+ Ax ^ +1 / 2 ~ S i-V2) 
U r = U* - ^ (F* +1 - FT) + ^ (D* +1/2 - D*_ 1/2 ) 

u " +1 = ^( U i + U D- 

Here D™ +1 / 2 is a simple finite difference approximation to D. 
Straightforward evaluation of S would be 



Sj+l/2 



( ° ^ 

\ Qj+1/2 + U j+l/2 S j+1/2 J 



(12) 



but we will see that some adjustment must be made. The approximation to the stochastic 
stress tensor, sj+1/2, is computed as 



4 + i/ 2 = ]J (v j+ iT j+ i + ViTj) H? +1/a (13) 

where V c is the volume of a cell and the K's are independent, Gaussian distributed random 
values with zero mean and unit variance. The approximation to the discretized stochastic 
heat flux, qj+1/2, is evaluated as 

ti+ip = y^ ( ^ +i(T ^ i)2+ ^ ( ^ )2) *w (14) 

These same stochastic flux approximations are used in all the continuum methods presented 
here. 

The stochastic components of the flux, Sj +1 / 2 , are independent, identically distributed 
Gaussian random variables with mean zero and variance a 2 for i — n, *. Substituting this 
into the MacCormack scheme we find that the variance in the flux at j + 1/2 is given by 



*A standard version of two-step Lax-Wendroff was also considered with similar but slightly poorer results. 
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( s G s * + H 2 ) = G) 2 ^( s ") 2 ) + G) 2<i(s * )2> 
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That is, the variance in the flux is reduced to half its original magnitude by the 
averaging used in the two-step MacCormack algorithm. We correct this effect by replacing 



2.2 Piecewise Parabolic Method 

In [37] a piecewise linear second-order Godunov scheme was shown to be effective for solv- 
ing the fluctuating Burgers' equation. We considered two versions of higher-order Godunov 
methods for the LLNS, a piecewise linear version [38] and the Piecewise Parabolic Method 
(PPM) introduced in [39]. The PPM algorithm, based on the direct Eulerian version pre- 
sented in [10], produced considerably better results than the piecewise linear scheme. Since 
our goal is to preserve fluctuations, we do not limit slopes and we do not include discontinuity 
detection in the algorithm. 

For this scheme the hyperbolic terms of the LLNS equations are considered in terms 
of hydrodynamic and local characteristic variables. In hydrodynamic variables we have 
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The local characteristic variables are interpolated via a fourth-order scheme to the 
left (— ) and right (+) edges of each cell: 

W£t = ^(L,V, + LjV j±1 ) - ^(L,V m + L,V i±2 ), (17) 

where Lj is the matrix whose rows are the left eigenvectors of A evaluated at Vj. 

These values, together with the cell-centered value W™ = LjV,-, are used to construct 
a parabolic profile Wj^(d) for each characteristic variable k in each cell, 

W(0) = W i( _ + OAWj + 6(1 - 9)W j6 , (18) 

where 



x 



(j - ±)Ax 



Ax 

AW" = W" + -W™_, and 
WJ 6 = 6(Wj-i(W^ + + W«J). 

Time-centered updates are based on the sign of each local characteristic wavespeed, 
W "+V2 = f i W,, fc (0) dB, ±A, fc > 



W™ ±fc otherwise 



where z/ i>fc = A i)fe ^. 

Finally, the time-centered values are transformed back into primitive variables and 
used as inputs to a Riemann problem at each cell edge. We use the approximate Riemann 
solver discussed in [41J. This approach iterates the phase space solution in the u — p plane, 
approximating the rarefaction curves by the Hugoniot locus. The overall approach is able 
to handle strong discontinuities and is second-order in wave strength. 

Approximations to the viscous and stochastic flux terms are discussed in section 12.11 
For our PPM algorithm we center the viscous update in time, so that the complete update 
is as follows: 

U* = un_^F» + ^(D» + Sy) (19) 

ur 1 = u " - if F " + \ ( {pi + s" + d; + s*) . (20) 

As discussed in section I2TT1 for the PPM scheme we use the adjusted stochastic flux approx- 
imation Sj = \/2Sj, since the averaging in the time-centering reduces the variance in the 
flux to half its original magnitude. 
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2.3 Variance-preserving third-order Runge-Kutta 

Equilibrium tests, presented in detail in the next section, show that neither stochastic version 
of the traditional numerical methods discussed above accurately represents the fluctuations 
in the LLNS equations. The principal difficulty arises because there is no stochastic forcing 
term in the mass conservation equation. Accurately capturing density fluctuations requires 
that the fluctuations be preserved in computing the mass flux. Another key observation is 
that the representation of fluctuations in the above schemes is also sensitive to the time step, 
with extremely small time steps leading to somewhat improved results. This suggests that 
temporal accuracy also plays a significant role in capturing fluctuations. Based on these 
observations we have developed a new discretization aimed specifically at capturing fluctua- 
tions in the LLNS equations. The method is based on a third order Runge-Kutta temporal 
integrator (RK3) combined with a centered discretization of hyperbolic and diffusive fluxes. 
The RK3 discretizaton can be written in the following three-stage form: 

u; +1/3 = u "-^(^Vi/ 2 --^-v 2 ) (2i) 

TT n+2/3 _ 3 1 n+1/3 _ 1 / At \ n+ i/ 3 n+ i/ 3 . 

u i ~ 4 U i + 4 u i A\AxJ [ j+1/2 i- 1 ' 2 ' 1 ' 

TT n+l _ l TT n j_ 2 TT n+2/3 2 / At \ n+2 /3 ^-n+2/3-, fOO\ 



where T = — F + D + S. 

Combining the three stages, we can write 



jjn+1 _ jjn 

3 ' 3 X,. 



( T n — T n \ _l_ ( •7T"+ 1 / 3 _ 7r«+l/3\ , z / -j^n+2/3 _ ^~n+2/3\ 
g^i+1/2 3-1/2) g^j+1/2 j-1/2 > + 3^+1/2 ^3-1/2 > 



The stochastic components of the flux, S™^ 2 are independent, identically distributed 
Gaussian random variables with mean zero and variance a 2 for t = 0, |, |. Substituting this 
into the combined update we find that the variance in the flux at j + 1/2 is given by 

I) 2 <(^° +1/2 ) 2 > + Q) 2 «< 3 1/2 ) 2 > + g) 2 (« 3 1/2 ) 2 ) 

O 2 

T" 

Thus, in the course of the RK3 algorithm, the variance in the flux is reduced to half 
its original magnitude, so again we replace Sj+1/2 by Sj+i/2 = v2S J - +1 / 2 , as discussed in 
section 12.11 and compute equations (12111231) using T = — F + D + S. 
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However, this treatment does not directly affect the fluctuations in density, since S 
does not appear in the continuity equation. We can correct this effect via a special interpo- 
lation scheme: by augmenting the variance to compensate for the density reduction arising 
from the temporal averaging, the fluctuations are preserved in the mass flux computation. 

We interpolate J (and the other conserved quantities) from cell-centered values: 

Jj+1/2 = ai{Jj + Jj+i) - a 2 (Jj-i + Jj+2), (24) 

where 

Ql = ( v / 7 + l)/4and (25) 

a 2 = (V7-l)/A. (26) 

Then in the case of constant J we have exactly Jj+1/2 = J and (SJj +1 , 2 ) = 2(5J 2 ), as 
desired; the interpolation is consistent and compensates for the variance-reducing effect of 
the multi-stage Runge-Kutta algorithm. The interpolation formula is similar to the PPM 
spatial construction except in the PPM construction ot\ = 7/12 and 02 = 1/12. Tests based 
on these alternative weights produced results intermediate to the RK3 scheme and the PPM 
scheme. We also considered interpolation of primitive variables but found that interpolation 
based on primitive variables led to stable but undamped oscillatory behavior. Finally, the 
diffusive terms D are discretized with standard second-order finite difference approximations. 

2.4 Boundary Conditions 

In sections |3] and H] we consider test problems for the various PDE algorithms on either 
a periodic computational domain, a computational domain bounded by thermal walls, or a 
computational domain bounded by infinite reservoirs. Boundary conditions are implemented 
using ghost cells. For the periodic and reservoir boundaries, it is straightforward to determine 
the ghost cell data. 

For the case of thermal walls, in addition to ghost cells we also use a one-sided finite 
difference formulation to approximate u x and T x in the calculation of the diffusive flux. The 
treatment of the hyperbolic flux at thermal walls varies by method. 

For thermal wall boundaries in MacCormack, conserved quantities are reflected across 
the boundaries of the domain. The temperature in the ghost cells is determined by linear 
extrapolation, and the no-flow condition is enforced by setting the velocity terms of the 
hyperbolic flux to zero within the ghost cells. 

For thermal wall boundaries in PPM, ghost cells are populated by reflecting primitive 
variable values across the domain boundaries, and the temperature in the ghost cells is deter- 
mined by linear extrapolation. The PPM routine takes as input the cell-centered primitive 
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variable data and returns a Riemann solution at each cell edge. On the domain boundaries, 
we modify these Riemann solutions by enforcing fixed wall temperature (i.e., the pressure 
at the wall is taken to be a function of the fixed wall temperature) before computing the 
hyperbolic flux across each edge. 

For thermal wall boundaries in RK3, conserved quantities are reflected across the 
boundaries of the domain and then interpolated onto cell edges. At the domain boundaries we 
employ a Riemann solver, which ensures that the boundary treatment respects characteristic 
compatibility relations at the physical boundaries. At the physical boundaries, the primitive 
variable values derived from the conserved-quantity interpolants are modified to enforce zero 
velocity and fixed wall temperature. This vector of primitive variables provides the input 
to the Riemann problem on the interior side of the boundary. The input to the Riemann 
problem on the exterior side of the boundary is the reflection of the interior input data. 
The treatment of reservoir boundaries is similar. However, ghost cells are populated with 
reservoir data, wall conditions are not enforced, and the input to the Riemann problem on 
the exterior side of the boundary is the reservoir data. 

3 Numerical Tests — Equilibrium 

This section presents results from a variety of scenarios in which the three schemes described 
above were tested. The physical domain is chosen to be compatible with DSMC particle 
simulations; see Table Q] for the system's parameters and Appendix B for a description of 
DSMC. The domain is partitioned into 40 cells of equal size Ax and hyperbolic and diffusive 
stability constraints determine the maximum time step At: 



where the sound speed c s = y 7-P/p, rj = T)(T), and k = k(T)] the overline indicates reference 
values (e.g., equilibrium values around which the system fluctuates). For the reference state 
(Argon at STP) and a cell width of Ax ~ 10~ 6 cm the time step used was At = 10 -12 s. 

3.1 Variances at equilibrium 

The first benchmark for our numerical schemes is recovering the correct variance of fluctu- 
ations for a system at equilibrium. For this initial test problem, we take a periodic domain 
with zero net flow and constant average density and temperature. Similar results, not pre- 
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System length 


1.25 x 10" 4 


Reference mean free path 


6.26 x 10" 6 


System volume 


1.96 x 10~ 16 


Time step 


1.0 x 10~ 12 


Number of cells 


40 


Number of samples 


10 7 


Number of DSM C particles 


5265 


DSMC collision grid size 


3.13 x 10" 6 



Table 1: System parameters (in cgs units) for simulations of a dilute gas in a periodic domain. 



sented here, were obtained for the case of constant non-zero net flow. The variances are 
computed in 40 spatial cells from 10 samples and then averaged over the cells. 

Table [2] compares the theoretical variances (see Appendix A) with those measured 
in the three stochastic PDE schemes and the DSMC particle simulation. The MacCormack 
and PPM schemes do relatively poor job (9 — 16% error) for the variances of density and 
energy. Better PPM results are obtained by decreasing our value of At by a factor of 10, 
to 10~ 13 . However, it is not desirable to run simulations at such a small time step. Only 
the third-order Runge-Kutta integrator generates the correct variance of density and energy 
while advancing with time steps near the stability limit. 



Table 2: Variance in conserved quantities at equilibrium (computed 
values are accurate to approximately 0.1%). 
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MacCormack scheme 
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13.31 


2.61 x 10 10 


Piece-wise Parabolic Method 
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13.27 


2.58 x 10 10 


Runge-Kutta (3 rd order) 


2.32 x 10- 


-8 


13.65 


2.87 x 10 10 


Molecular simulation (DSMC) 


2.35 x 10- 


-8 


13.21 


2.79 x 10 10 


Percentage difference (MacCormack) 


-14.3% 




2.3% 


-9.3% 


Percentage difference (PPM) 


-16.0% 




2.0% 


-10.3% 


Percentage difference (RK3) 


-1.3% 




4.9% 


-0.1% 


Percentage difference (DSMC) 


0.0% 




1.6% 


-3.1% 
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3.2 Spatial correlations at equilibrium 

Figures [TH3] depict the spatial correlation of conserved variables, that is, (SpjSpj*), (5Jj5Jj*), 
and (5Ej5Ej*), where j* is located at the center of the domain. These figures show results 
computed by the MacCormack, PPM, and RK3 schemes, along with the theoretical values 
of the correlations (see Appendix A) and molecular simulation data (see Appendix B). For 
the MacCormack and PPM schemes the spatial correlations of density fluctuations and 
energy fluctuations have significant spurious oscillations near the correlation point (see Figs.[T] 
and [3]). All three schemes do well in reproducing the expected correlations of momentum 
fluctuations. Figure 0] depicts (5pj5Jj*), which has a theoretical value of zero since the net 
flow is zero; all three schemes correctly reproduce this result. 



x Dsyc 
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x/L -2 1 ' ' ' ' ' ' ' ' ' ' 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x/L 



Figure 1: Spatial correlation of density fluc- 
tuations. Solid line is(5pi5pj) = {5p 2 )8^ (see 
equations (13811391)). 



Figure 2: Spatial correlation of momentum 
fluctuations. Solid line is(5Ji5Jj) = {5J 2 )5^ 
(see equations ( 1MI471) ). 



3.3 Time correlations at equilibrium 

The time correlation of density fluctuations is of interest because its temporal Fourier trans- 
form gives the spectral density, which is measured experimentally from light scattering spec- 
tra [121 H3] . From the LLNS equations, this time correlation can be written as 

(5p(w,t)5p(w,t + r)) / 1\ 2 1 2 

TTTT \w = 1 ~ - exp{~w 2 D T r} + - exp{-w 2 Tr} cos(c s wr) 

ST — D 

H -wexp{— w 2 Tt} sin(c s wr) (29) 

7 c s 
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Figure 3: Spatial correlation of energy fluctu- 
ations. Solid line ^(SEiSEj) = (5E 2 }5^ (see 
equations (HBl HSjl ). 



Figure 4: Spatial correlation of density- 
momentum fluctuations. 



where w = 2irn/L is the wavenumber, 7 = c p /c v is the ratio of specific heats, = n/pc v 
is the thermal diffusivity, D v = -^fj/p is the longitudinal kinematic viscosity, c s is the sound 
speed, and T = \[D V + (7 — 1)Dt] is the sound attenuation coefficient. 

In our numerical calculations the density is represented by cell averages p iy i = 1, . . . , M c , 
and the time correlation is estimated from the mean of N samples, 

1 - 

(6p(w,t)5p(w,t + r)) N =- R(t)R(t + r) (30) 



N 

samples 



with 

M, 



R(t) = ^ r J2p i M^nxi/L). (31) 

We have 



c i=i 



(5p(w,t)6p(w,t + t)} = lim (6p(w, t)6p(w, t + t))n- (32) 
From the above we find the normalization of the theoretical result may be expressed as 

(5p 2 (w, £)} = (R(t) 2 ) = — ^(SpiSpj) m^wnxi/L) smfinnxj/L) 

c i=i j=i 

= M. (33) 
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We restrict our attention to the lowest wavenumber (i.e., n = 1) because for the system 
sizes we consider the theoretical result, (1291) . is not accurate at short wavelengths due to 
mean-free-path corrections. 

In the left-hand panel of figure we present time correlation results from our equi- 
librium problem on a periodic domain. We compare results from the MacCormack, PPM, 
and RK3 methods with the theoretical time correlation, equation (129]) . and with molecular 
simulation data (see Appendix A). We find reasonable agreement among all the results, up to 
the time when a sound wave has crossed the system (k4x 10~ 9 seconds). Due to finite size 
effects the theory is only accurate for short times but the agreement among the numerical 
PDE schemes and DSMC molecular simulation is good. 

The right-hand panel of figure [5] shows time correlation results for the equilibrium 
problem on a domain with thermal walls rather than periodic boundaries; we find good 
agreement for this problem as well, at least for times less than the sound crossing time. 
For later times, the time correlation is sensitive to the acoustic impedance of the thermal 
wall. For this case, MacCormack under-predicts the correlation at early time while PPM 
shows significant deviation near t = 5 x 1CT 8 . Both MacCormack and the RK3 scheme 
deviate somewhat from DSMC at late time. Overall, however, the RK3 scheme captures the 
temporal correlation better than either of the other two PDE schemes. 

4 Numerical Tests — Non-equilibrium 

The results from the section above indicate that of the three stochastic PDE schemes, the 
third-order Runge-Kutta method (RK3) consistently out-performs the other two schemes. 
In this section we consider two more numerical tests, spatial correlations in a temperature 
gradient and diffusion of a standing shock wave, but restrict our attention to the RK3 scheme, 
comparing it with DSMC molecular simulations. 

4.1 Spatial correlations in a temperature gradient 

In the early 1980 's, a variety of statistical mechanics calculations predicted that a fluid 
under a non-equilibrium constraint, such as a temperature gradient, would exhibit long- 
range correlations of fluctuations [H]. Furthermore, quantities that are independent at 
equilibrium, such as density and momentum fluctuations, also have long-ranged correlations. 
These predictions were qualitatively confirmed by light scattering experiments [45j, yet the 
effects are subtle and difficult to measure accurately in the laboratory. Molecular simulations 
confirm the predicted correlations of non-equilibrium fluctuations for a fluid subjected to a 
temperature gradient [16], [23] and to a shear [17] . 
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Figure 5: Time correlation of density fluctuations for equilibrium prob- 
lem, on a periodic domain (left panel) and a domain with specular wall 
boundaries (right panel). 



We consider a system similar to that of section [3T31 but with a temperature gradient. 
Specifically, the boundary conditions are thermal walls at 273K and 819K. Figure [6] shows 
the correlation of density and momentum fluctuations measured in an RK3 calculation and 
by DSMC simulations. The two sets of data are in good agreement and are in agreement 
with earlier work on this problem [16J [23]. The major discrepancy is the under-prediction 
of the negative peak correlation near j*. Extensive tests suggest that this effect is hard 
to capture with a continuum solver because of the tension between variance reduction and 
spatial correlations in computing the mass flux at cell edges from cell-centered data. 

4.2 Random Walk of a Standing Shock 

In our final numerical study we consider the random walk of a standing shock wave due to 
spontaneous fluctuations. Shock diffusion is well-known in other particle simulations, such 
as shock tube modeling by DSMC, which must correct for the drift when measuring profiles 
for steady shocks. [H] The general problem has been also been analyzed for simple lattice 
gas models (33 [501 ED 1521 EZ! ■ 
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Figure 6: Spatial correlation of density and mo- 
mentum fluctuations for a system subjected to 
a temperature gradient. Compare with Fig. HI 



System length 
System volume 
Number of cells 


5 x 1(T 4 
7.84 x 1(T 16 
160 


Reference mean free path 
Time step 
Mach number 


6.26 x 10~ 6 
1.0 x 10~ 12 
2.0 


RHS mass density 
RHS velocity 
RHS temperature 
RHS sound speed 


1.78 x 1(T 3 
-61562 

273 
30781 


LHS mass density 
LHS velocity 
LHS temperature 
LHS sound speed 


4.07 x 10~ 3 
-26933 

567 
44373 



Table 3: System parameters (in cgs units) for simulations of a standing shock, Mach 2.0 



Mass density and temperature on the right-hand side of the shock are given the same 
values as in our equilibrium problem; values of density and temperature on the left-hand side 
are derived from the Rankine-Hugoniot relations. The velocity on both sides of the shock 
are specified to satisfy the Rankine-Hugoniot conditions and to make the unperturbed shock 
wave stationary in the computational domain. We consider three different shock strengths, 
Mach 2, Mach 1.4, and Mach 1.2 (see table [3]). The boundary treatment consists of infinite 
reservoirs with the same states as the initial conditions. For this test problem we use a 
longer computational domain, in order to capture (unlikely) shock drift of several standard 
deviations. 

Here we focus on the variance of the shock location as a function of time. We define 
a shock location for density, u p (t) by fitting a Heaviside function to the integrated density, 
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o(t) pL/2 t-L/2 

pidx+ PrcLx = I p(x,t)dx . (34) 

•Zr/2 Jcr p (t) J-L/2 

Solving for <r p (t) gives 

cr p (t) = L s (35) 

Pl- Pr 

L /2 

where p = L~ x §_ L , 2 p(x, t) dx is the instantaneous average density. The shock location for 
pressure, o>, is analogously defined. We estimate & p (t) and cr p (t) as functions of time from 
ensembles of 4000 simulations. For the PDE simulations, we initialize with discontinuous 
shock profiles. One would expect the shock location to fluctuate with a diffusion similar to 
that of a simple random walk [51 J, so averaging over ensembles from the same initial state 
we would expect to find 

(5a 2 p ) « 2 V p t and (6o%) « 2 V p t (36) 

with shock diffusion coefficients, T> p and T> P) that depend on shock strength. Note that this 
expression for the variance is not accurate at very short times (due to transient relaxation 
from the initial state) or at very long times (due to finite system size). 

Figure [7] shows results for the variance in the shock position from an ensemble of runs 
versus time. After the initial transients, the slopes are constant with the strongest shocks 
exhibiting the least drift (T> = (Ma — 1) _1 ) and with a p and o> giving similar diffusion 
coefficients. DSMC data is initially noisy so it has different initial transients and "diffuses" 
farther than the PDE. However, after the transients, the DSMC and the RK3 simulations 
have essentially the same slope, as a function of Mach number. This indicates that the 
third-order Runge-Kutta scheme is accurately capturing the shock-drift random walk. 



5 Summary and Concluding Remarks 

In this paper we develop and analyze several finite- volume schemes for solving the fluctuat- 
ing Landau-Lifshitz compressible Navier-Stokes equations in one spatial dimension. Methods 
based on standard CFD discretizations were found not to accurately represent fluctuations in 
an equilibrium flow. We have introduced a centered scheme based on interpolation schemes 
designed to preserve fluctuations combined with a third-order Runge-Kutta (RK3) tem- 
poral integrator that was able to capture the equilibrium fluctuations. Further tests for 
non-equilibrium systems confirm that the RK3 scheme correctly reproduces long-ranged cor- 
relations of fluctuations and stochastic drift of shock waves, as verified by comparison with 
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Figure 7: Variance of shock location for mass density profile (left panel, (5a p (t) 2 )) and 
pressure profile (right panel, (5<jp(t) 2 )). Estimated variances (4000-run ensembles) versus 
time t for a deterministically steady shock of Mach number 1.2, 1.4, or 2.0. Solid lines are 
for RK3, dashed lines are from DSMC molecular simulations. 

molecular simulations. It is worth emphasizing that the ability of continuum methods to ac- 
curately capture fluctuations is fairly sensitive to the construction of the numerical scheme. 
Minor variations in the numerics can lead to significant changes in stability, accuracy, and 
behavior. 

The work discussed here suggests a number of additional studies. Further analysis 
is needed on the treatment of thermal and reservoir boundary conditions. The methods 
here can also be extended to three dimensions (for which the stochastic stress tensor is 
more complex) and we can include concentration as a hydrodynamic variable to allow the 
methodology to be applied to a number of other flow problems. Finally, we are integrating 
our new stochastic PDE solver into our existing Adaptive Mesh and Algorithm Refinement 
(AMAR) programs [53J. A stochastic AMAR simulation will not only model hydrodynamic 
fluctuations at multiple grid scales but will, by incorporating DSMC simulations at the finest 
level of algorithm refinement, also capture molecular- level physics. 
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Appendix A: Equilibrium Fluctuations 



For infinite systems, at thermodynamic equilibrium both conserved and hydrodynamic vari- 
ables are spatially uncorrelated at equal times. For example, 

(6pt(t)6 Pj (t)) = (5p 2 )5«. (37) 

For conserved variables there is a finite size correction, specifically, 

(6 Pt (t)6p 3 (t)} = (i - iic W>*5 - acWXi - *§) (38) 

for z,j = 1,...,M C , where M c is the number of cells in the system. The variances are 
well-known from equilibrium statistical mechanics (§112, [54]). 

The variance of mass density depends on the compressibility (i.e., the equation of 
state) of the fluid. In general, 

W)=f^ (39) 

where N c and (SN 2 ) are the mean and variance of the number of particles in a cell. We 
calculate N c = ~pV c / m, where V c is the volume of a cell and m is the mass of a particle. For 
an ideal gas N c is Poisson distributed so (SN 2 ) = N c and (Sp 2 ) = p 2 /N c . The more general 
result is (SN 2 ) = axphsT N c /m where «t is the isothermal compressibility. 
The variances of fluid velocity and temperature in a cell are 

(gy?) = ^ = 2 and (40) 
Wc N c 



c v pV c c v N c 



{6T 2 ) = M1 = -^E, (41) 



where Ct = yksT/m is the thermal speed (and the standard deviation of the Maxwell- 
Boltzmann distribution). The covariances are (5pSu) = (5p5T) = (5u5T) = 0. 

The variances and covariances of the mechanical densities at equilibrium are 

{5p5J} = pJA p (42) 

(6p6E) = pEA p (43) 

(8 J 2 ) = J 2 A p + p 2 C 2 A u (44) 

(5J5E) = JEA P + J pClA u (45) 

(5E 2 ) = E 2 A p + J 2 C 2 A u + clp 2 T 2 A T (46) 

where A p = (8p 2 }/p 2 , A u = (5u 2 )/C%, and A T = (5T 2 }/T 2 . For a dilute gas A p = A u = 
1/N C , and = 2/(3N c ). Again, corrections must be made for conserved quantities in the 
case of a finite domain: 

(SJiWJtf)) = (1 - M~ 1 )(5J 2 )5-j - M- 1 (5J 2 )(1 - tfg), (47) 
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(SEkWW)) = (1 - M^)(6E 2 )5« - M-i(5E*)(l - (48) 

Appendix B: DSMC Simulations 

The algorithms presented here for the stochastic LLNS equations were validated by com- 
parison with molecular simulations. Specifically, we used the direct simulation Monte Carlo 
(DSMC) algorithm, a well-known method for computing gas dynamics at the molecular scale; 
see [55j EH] f° r pedagogical expositions on DSMC, [H] for a complete reference, and [57] for 
a proof of the method's equivalence to the Boltzmann equation. As in molecular dynamics, 
the state of the system in DSMC is given by the positions and velocities of particles. In each 
time step, the particles are first moved as if they did not interact with each other. After 
moving the particles and imposing any boundary conditions, collisions are evaluated by a 
stochastic process, conserving momentum and energy and selecting the post-collision angles 
from their kinetic theory distributions. DSMC is a stochastic algorithm but the statistical 
variation of the physical quantities has nothing to do with the "Monte Carlo" portion of the 
method. The equilibrium and non-equilibrium variations in DSMC are the physical spectra 
of spontaneous thermal fluctuations, as confirmed by excellent agreement with fluctuating 
hydrodynamic theory [3H [23] and molecular dynamics simulations [58l 124]. 

The simulated physical system is a dilute monatomic hard-sphere gas in a rectangular 
volume with periodic boundary conditions in the y and z directions. The boundary conditions 
in the x direction are either periodic, specular (i.e., elastic reflection of particles), or a pair 
of parallel thermal walls. The physical parameters used are presented in Table HI Samples 
are taken in forty rectangular cells perpendicular to the x-direction. 

References 

[1] R.K. Pathria. Statistical Mechanics. Butterworth-Heinemann, Oxford, 1996. 

[2] G. Karniadakis, A. Beskok, and N. Aluru. Microflows and Nanoflows : Fundamentals 
and Simulation. Springer, New York, 2005. 

[3] CM. Ho ; Y.C. Tai. Micro-electro-mechanical systems (MEMS) and fluid flows. Annu. 
Rev. Fluid Mech., 30:579-612, 1998. 

[4] M. Gad el Hak. The fluid mechanics of microdevices -the Freeman Scholar lecture. J. 
Fluids Eng., 121:5-33, 1999. 

[5] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter. Molecular Biology 
of the Cell. Garland, New York, 4th edition, 2002. 



21 



[6] R. D. Astumian ; P. Hanggi. Brownian motors. Physics Today, pages 33-39, November 
2002. 

[7] G. Oster. Darwin's motors. Nature, 417:25, 2002. 

[8] R.K. Soong; G.D. Bachand; H.P. Neves; A.G. Olkhovets; H.G. Craighead; CD. Monte- 
magno. Powering an inorganic nanodevice with a biomolecular motor. Science, 290:1555, 
2000. 

[9] T.Y. Tsong. Na,K-ATPase as a Brownian motor: Electric field-induced conformational 
fluctuation leads to uphill pumping of cation in the absence of ATP. J. Bio. Phys., 
28:309-325, 2002. 

[10] H.G. Craighead. Nanoelectromechanical systems. Science, 290:1532, 2000. 

[11] A. van Oudenaarden and S.G. Boxer. Brownian ratchets: Molecular separations in lipid 
bilayers supported on patterned arrays. Science, 285:1046-48, 1999. 

[12] J. Bader, R. Hammond, S. Henck, M. Deem, G. McDermott, J. Bustillo, J. Simpson, 
G. Mulhern, and J. Rothberg. DNA transport by a micromachined Brownian ratchet 
device. Proc. Natl. Acad. Set., 96:13165-9, 1999. 

[13] C. Van den Broeck, R. Kawai, and P. Meurs. Exorcising a Maxwell demon. Phys. Rev. 
Lett, 93:090601, 2004. 

[14] P. Meurs, C. Van den Broeck, and A.L. Garcia. Rectification of thermal fluctuations in 
ideal gases. Phys. Rev. E, 70:051109, 2004. 

[15] B. Nowakowski and A. Lemarchand. Sensitivity of explosion to departure from partial 
equilibrium. Physical Review E, 68:031105, 2003. 

[16] A. Lemarchand and B. Nowakowski. Fluctuation-induced and nonequilibrium-induced 
bifurcations in a thermochemical system. Molecular Simulation, 30(1 1-12) :773-780, 
2004. 

[17] L.D. Landau and E.M. Lifshitz. Fluid Mechanics, volume 6 of Course of Theoretical 
Physics. Pergamon, 1959. 

[18] M. Bixon and R. Zwanzig. Boltzmann-Langevin equation and hydrodynamic fluctua- 
tions. Phys. Rev., 187(l):267-272, Nov 1969. 

[19] R. F. Fox and G. E. Uhlenbeck. Contributions to non-equilibrium thermodynamics. I. 
Theory of hydrodynamical fluctuations. Phys. Fluids, 13(8):1893-1902, 1970. 



22 



[20] G.E. Kelly and M.B. Lewis. Hydrodynamic fluctuations. Physics of Fluids, 14(9):1925- 
1931, 1971. 

[21] E. Calzetta. Relativistic fluctuating hydrodynamics. Class. Quantum Grav., 15:653, 
1998. 

[22] P. Espa nol. Stochastic differential equations for non-linear hydrodynamics. Physica A, 
248:77, 1998. 

[23] M. Malek-Mansour, A.L. Garcia, G.C. Lie, and E. Clementi. Fluctuating hydrodynamics 
in a dilute gas. Phys. Rev. Lett, 58:874-877, 1987. 

[24] M. Mareschal, M. Malek-Mansour, G. Sonnino, and E. Kestemont. Dynamic struc- 
ture factor in a nonequilibrium fluid: A molecular-dynamics approach. Phys. Rev. A, 
45:7180-7183, May 1992. 

[25] A.J.C. Ladd. Short-time motion of colloidal particles: Numerical simulation via a 
fluctuating lattice-Boltzmann equation. Phys. Rev. Lett, 70(9): 1339-1342, Mar 1993. 

[26] N. Sharma and N.A. Patankar. Direct numerical simulation of the Brownian motion of 
particles by using fluctuating hydrodynamic equations. J. Comput. Phys., 201(2):466- 
486, 2004. 

[27] M. Moseler and U. Landman. Formation, stability, and breakup of nanojets. Science, 
289(5482):1165-1169, 2000. 

[28] H.P. Breuer and F. Petruccione. A master equation description of fluctuating hydrody- 
namics. Physica A, 192:569-588, February 1993. 

[29] H.P. Breuer and F. Petruccione. A master equation approach to fluctuating hydrody- 
namics: Heat conduction. Phys. Lett. A, 185:385-389, February 1994. 

[30] M. Serrano and P. Espa nol. Thermodynamically consistent mesoscopic fluid particle 
model. Phys. Rev. E, 64(4):046115, Sep 2001. 

[31] M. Grmela and H.C. Ottinger. Dynamics and thermodynamics of complex fluids. I. 
Development of a general formalism. Phys. Rev. E, 56(6):6620-6632, Dec 1997. 

[32] G. De Fabritiis, P.V. Coveney, and E.G. Flekky. Multiscale dissipative particle dynam- 
ics. Philos. Trans. R. Soc. London, Ser. A, 360:317-331, 2002. 

[33] M. Serrano, G. De Fabritiis, P. Espa nol, E.G. Flekk0y, and P.V. Coveney. Mesoscopic 
dynamics of Voronoi fluid particles. J. Phys. A, 35(7):1605-1625, 2002. 



23 



[34] A.L. Garcia, M. Malek-Mansour, G. Lie, and E. Clementi. Numerical integration of the 
fluctuating hydrodynamic equations. J. Stat. Phys., 47:209, 1987. 

[35] F.J. Alexander, A.L. Garcia, and D.M. Tartakovsky. Algorithm refinement for stochastic 
partial differential equations: I. Linear diffusion. J. Comput. Phys., 182(l):47-66, 2002. 

[36] F.J. Alexander, A.L. Garcia, and D.M. Tartakovsky. Algorithm refinement for stochastic 
partial differential equations: II. Correlated systems. J. of Comp. Phys., 207:769-787, 
2005. 

[37] J.B. Bell, J. Foo, and A. Garcia. Algorithm refinement for the stochastic Burgers' 
equation. J. Comp. Phys., (in press), 2006. 

[38] P. Colella. A direct Eulerian MUSCL scheme for gas dynamics. SIAM J. Sci. Stat. 
Comput, 6:104-117, 1985. 

[39] P. Colella and P.R. Woodward. The Piecewise Parabolic Method (PPM) for gas- 
dynamical simulations. J. of Comp. Phys., 54:174-201, 1984. 

[40] R.E. Miller and E.B. Tadmor. The quasicontinuum method: overview, applications and 
current directions. J of Comput. Aided Mater. Des., 9 (3): 203-39, 2002. 

[41] P. Colella and H.M. Glaz. Efficient solution algorithms for the Riemann problem for 
real gases. J. of Comp. Phys., 59:264-289, 1985. 

[42] B. J. Berne and R. Pecora. Dynamic Light Scattering: With Applications to Chemistry, 
Biology, and Physics. Dover, 2000. 

[43] J. P. Boon and S. Yip. Molecular Hydrodynamics. Dover, 1991. 

[44] R. Schmitz. Fluctuations in nonequilibrium fluids. Physics Reports, 171:1, 1988. 

[45] D. Beysens, Y. Garrabos, and G. Zalczer. Experimental evidence for Brillouin asym- 
metry induced by a temperature gradient. Phys. Rev. Lett., 45:403, 1980. 

[46] A.L. Garcia. Nonequilibrium fluctuations studied by a rarefied gas simulation. Phys. 
Rev. A, 34:1454, 1986. 

[47] A.L. Garcia, M. Malek-Mansour, G.C. Lie, M. Mareschal, and E. Clementi. Hydrody- 
namic fluctuations in a dilute gas under shear. Phys. Rev. A, 36:4348-4355, 1987. 

[48] G.A. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon, 
Oxford, 1994. 



24 



[49] F.J. Alexander, Z. Cheng, S.A. Janowsky, and J.L. Lebowitz. Shock fluctuations in the 
two-dimensional asymmetric simple exclusion process. J. Stat. Phys., 68(5-6):761-785, 
1992. 

[50] F.J. Alexander, S.A. Janowsky, J.L. Lebowitz, and H. van Beijeren. Shock fluctuations 
in one-dimensional lattice fluids. Phys. Rev. E, 47:403-410, 1993. 

[51] P.A. Ferrari and L.R.G. Fontes. Shock fluctuations in the asymmetric simple exclusion 
process. Probab. Theory Related Fields, 99(2):205-319, 1994. 

[52] S.A. Janowsky and J.L. Lebowitz. Finite-size effects and shock fluctuations in the 
asymmetric simple-exclusion process. Phys. Rev. A, 45:618-625, January 1992. 

[53] A.L. Garcia, J.B. Bell, W.Y. Crutchfield, and B.J. Alder. Adaptive mesh and algorithm 
refinement using Direct Simulation Monte Carlo. J. Comput. Phys., 154(1):134-155, 
1999. 

[54] L.D. Landau and E.M. Lifshitz. Statistical Physics, volume 5 of Course of Theoretical 
Physics. Pergamon, third ed., part 1 edition, 1980. 

[55] F.J. Alexander and A.L. Garcia. The Direct Simulation Monte Carlo method. Com- 
puters in Physics, ll(6):588-593, 1997. 

[56] A.L. Garcia. Numerical Methods for Physics. Prentice Hall, 2nd edition, 2000. 

[57] W. Wagner. A convergence proof for Bird's Direct Simulation Monte Carlo method for 
the Boltzmann equation. J. Stat. Phys., 66:1011, 1992. 

[58] M. Malek-Mansour, A.L. Garcia, J.W. Turner, and M. Mareschal. On the scattering 
function of simple fluids in finite systems. J. Stat. Phys., 52:295, 1988. 



25 



